subroutine rgp_fortran_energy_aaaa( &
       no, & ! (in)    number of occupied orbitals
       nv, & ! (in)    number of unoccupied orbitals
      Tia, & ! (in)    one electron integrals
    Tiajb, & ! (in)    two electron integrals in 1122 order
      Bkl, & ! (in)    inverse matrix
      Via, & ! (in)    new rows
     VBia, & ! (out)   contraction of B and V
    TVBjb, & ! (out)   contraction of T, V, and B
    TVBaj, & ! (out)   contraction of T, V, and B
   retval)   ! (inout) running total for the local energy

implicit none

integer,      intent(in)    ::     no
integer,      intent(in)    ::     nv
real(kind=8), intent(in)    ::    Tia(no,nv)
real(kind=8), intent(in)    ::  Tiajb(no,nv,no,nv)
real(kind=8), intent(in)    ::    Bkl(no,no)
real(kind=8), intent(in)    ::    Via(no,nv) ! second index labels which row
real(kind=8), intent(out)   ::   VBia(no,nv) ! second index labels which row
real(kind=8), intent(out)   ::  TVBjb(no,nv)
real(kind=8), intent(out)   ::  TVBaj(nv,no)
real(kind=8), intent(inout) :: retval

integer :: a, j

!real(kind=8) :: tempE

include 'formic/fqmc/rgp_fortran_interface.fpp'

! compute  VBia
call rgp_fortran_build_vbia(no, nv, Bkl, Via, VBia)

! compute TVBjb
call rgp_fortran_build_tvbjb(no, nv, Tiajb, VBia, TVBjb)

! compute TVBaj
call rgp_fortran_build_tvbaj_aaaa(no, nv, Tiajb, VBia, TVBaj)

! one electron integral term
retval = retval + ddot(no*nv, Tia, 1, VBia, 1)

!tempE = 0.00d+00

! add two electron "diagonal" term
retval = retval + 0.50d+00 * ddot(no*nv, TVBjb, 1, VBia, 1)
!tempE = tempE   + 0.50d+00 * ddot(no*nv, TVBjb, 1, VBia, 1)

! add two electron "off-diagonal" term
do j = 1,no
do a = 1,nv
  retval = retval - 0.50d+00 * TVBaj(a,j) * VBia(j,a)
  !tempE = tempE - 0.50d+00 * TVBaj(a,j) * VBia(j,a)
enddo
enddo

!retval = retval + tempE

!write(*, "(F30.12)") tempE
!call flush()

end subroutine rgp_fortran_energy_aaaa
